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Abstract. A generalized Black-Scholes equation is considered on the 
semi-axis. It is transformed on the interval (0, 1) in order to make the 
computational domain finite. The new parabolic operator degenerates 
at the both ends of the interval and we are forced to use the Garding 
inequality rather than the classical coercivity. A fitted finite volume el- 
ement space approximation is constructed. It is proved that the time 6- 
weighted full discretization is uniquely solvable and positivity-preserving. 
Numerical experiments, performed to illustrate the usefulness of the 
method, are presented. 
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1 Introduction 

The famous equation, proposed by F. Black, M. Scholes and R. Merton, see 
[13115] . is 

f + ^ 2 S^(g-^)-0,^(0,ooM e[ 0,T). 

This is a typical example of a degenerate parabolic equation [14]. It is well 
known [13119] that it can be transformed to the heat equation that allows us 
to overcome the degeneracy at S — 0. Many numerical methods, based on clas- 
sical finite difference schemes, applied to constant coefficients heat equations, 
are developed [1115] , However, when the problem has space-dependant coeffi- 
cients <7 and r one can not transform the Black-Scholes equation to a standard 
heat equation. Finite difference and finite element methods have been applied 
in |2I3I4I5I6I7I8I10I16I21| in order to solve this type of generalized Black-Scholes 
equations. In [TT] cubic B-splines are implemented. Often, the convergence of 
the full discretization is verified by numerical examples only. 

An effective method, that resolves the degeneracy, is proposed by S. Wang [20 
for the Black-Scholes equation with Dirichlet boundary conditions. The method 



is based on a finite volume formulation of the problem, coupled with a fitted 
local approximation to the solution and an implicit time-stepping technique. 
The local approximation is determined by a set of two-point boundary value 
problems (BVPs), defined on the element edges. This fitted technique originates 
from one-dimensional computational fluid dynamics [T5] . 

A modification of the discretization, originally presented in [20 , was proposed 
by L. Angermann |2] such that the method adequately treats the proper (natural) 
boundary condition at x = 0. Similar space discretization is derived in j6] for a 
degenerate parabolic equation in the zero-coupon bond pricing. 

The domain of S is the half real line. For numerical computation it is desirable 
to have a finite computational domain. The transformation in the next section 
converts S £ (0, oo) to x £ (0,1), decreasing significantly the computational 
costs. Also, for a call option, the solution V(S,t) is not bounded and from the 
numerical methods' point of view the problem transformation on a finite interval 
is better. The resulting equation has variable coefficients but this is not an 
essential difficulty for the numerical computation. However, the transformed 
equation degenerates at both ends of the finite interval. 

The present paper deals with a degenerate parabolic equation, Q, derived 
after transformation of the generalized Black-Scholes equation |l| to a finite 
interval. The degeneration at the both ends of the interval does not allow the use 
the Poincare-Friedrichs inequality and we are forced to investigate the differential 
problem with the Garding inequality rather than classical coercivity [5]. 

This paper is organized as follows. The model problem is presented in Section 
2, where we discuss our basic assumptions and some properties of the solution. 
The space discretization method is developed in Section 3. Section 4 is devoted to 
the time discretization, where we show that the system matrix on each time-level 
is an M-matrix so that the discretization is monotone. Numerical experiments 
are discussed in the last section. 

Some results, concerning the case of the transformed Black-Scholes equation 
above, are reported in |17| . 



2 The transformed problem 

We consider the generalized Black-Scholes equation |13I19| : 



~m + ^^ds^ + i^s-Dis^t))—- r(t)v = o, (s,t) e (o, oo) x (o,t), 

(1) 

where a — o~(t) denotes the volatility of the asset, r = r(t) is the risk-free 
interest rate, D = D(S,t) denotes the dividend of the dividend-paying asset. We 
also introduce d = d(S,t) such that D(S,t) = Sd(S,t) 1 where the dividend rate 
d = d(s,t) is continuously differentiable with respect to S. 

There are various choices for the final (payoff) condition, depending on the 
models. In the case of vanilla European option 



, „ , _ J max(S — E), for a call option, , , 

^ ' ' 1 max(£ — S 1 ), for a put option, ^ ' 

where E is the strike price. Another example is the bullish vertical spread payoff, 
defined by 

V(S,T)=max(S-E 1 )-max(S-E 2 ), (3) 

where E\ and E^ are two exercise prices, satisfying < E\ < E^. This represents 
a portfolio of buying one call option with exercise price E\ and issuing one call 
option with the same expiry date but a larger exercise price E^. For detailed 
discussion on this, we refer to |19| . 

We introduce the transformation [3T] 



s + P m s + P n 

The constant P rn is called mesh parameter. It controls the distribution of the 
mesh nodes w.r.t. S on the interval (0, oo). The higher the value of S, that 
we are interested in, the higher value of P m should be in order to obtain a 
reasonable accuracy. In the case of a call option, because of the nature of the 
terminal condition, P rn should be equal to E. 

The inverse transformation is S = P m x/(1 — x) and after plugging it in the 
Black-Scholes equation, (JT|) we obtain: 

+ ((1 - x)r{t) + xd{x, t))u = 0, {x, t)en T = Ox(0,T),O = (0, 1). 

We return to the original notation of the variable t for the sake of simplicity. 
The initial data for a call option reads 

u(x, 0) = uq{x) = max(2x — 1,0). (5) 

Being different from the classical parabolic equations, in which the principle 
coefficient is assumed to be strictly positive, the parabolic equation Q belongs 
to the second order differential equations with non-negative characteristic form 
1 14] . The main difficulty of such kind of equations is the degeneracy |18| . It can 
be easily seen that at x — and 1 = 1 ffl degenerates to 



du 



. . .„ . du 
= -r(t)u(0,t), — 



-d(l,t)u(l,t). 



It is well known by the Fichera theory for degenerate parabolic equations 
[T3] that at the degenerate boundaries x = and x = 1 the boundary conditions 
should not be given. 

For theoretical analysis of our discrete problem as well as for the construction 
of a fitted finite volume mass-lumping discretization we need to consider weak 
solutions of p}. We shall use the standard notations for the function space 



C m (fi) and C m (f2) of which a function and it's derivatives up to order m are 
continuous on Q (respectively J?). The space of square-integrable functions we 
denote by L 2 (Q) with the usual L 2 -norm || • || and the inner product (•, •). We also 
use the function space L°°(fi) with the norm || • W^. To handle the degeneracy 
in Q , we introduce the following weighted L^,-norm 

\\v\\ 0<w :=([ x 2 (l-x) 2 v 2 dx)^ 2 



Jo 

with corresponding inner product (u,v) w . Using L 2 (J?) and L^,(i2), we define 
the weighted Sobolev space H^(f2) := {v : v <E L 2 (Q), v' £ where v' 

denotes the weak derivative of v. Let || ■ ||x,«, be the functional on H^Q), defined 
by |M|i,u> = (IMIo + ll u 'llo.tu) 1/ ' 2 - Then it is easy to see that || • \\i^ w is a norm on 
H 1 ^); it is called weighted i/^-norm on H^fi). Furthermore, using the inner 
products in L 2 (Q) and L^(i?), we define a weighted inner product on H^(f2) by 
(•, -)h ■— (•, •) + (•, -) w and, consequently, the pair (•, •)#-) is a Hilbert 

space. Also, H^(f2) contains the conventional Sobolev space as a proper 

subspace. 

We rewrite Q in divergent form 

^ - - x) (a(x,t)®^ + b(x,t)i?j^J + c(x,t)u = 0, 

a(x, t) = ^<r 2 (t)x(l - x), b(x, t) = (r(t) - d(x, t) + a 2 (t)(2x - 1)) , (6) 

dd 

c(x,t) = (2-3x)r(t) - (6x 2 - 6x + l)a 2 (t) - (1 - 3x)d(x,t) - x(l - x) — . 

ox 

Let us introduce for w,v £ H^(f2) the bilinear form 

A(w, v; t) := (aw' + bw, v') + (cw, v) = (x(l — x)p(w), v) + (cw, v). 

Here the notation w' — ^ is used and the function p(w) — aw' + bw is 
the weighted flux density associated with w. We are in position to state the 
variational formulation of Q: 

Find u(t) € H^Q), such that for all v G H^(Q) 

^^-,vj+A(u(t),v;t)=Oa..e.m(0,T)axid (;0)=u . (7) 

The following result provides the weak coercivity and continuity of the bilinear 
form A(-, -,t). 

Lemma 1. Let w,v e C 1 ([0, 1]). Then: 
1. there exist constants C\ > and Ci > such that 

\A(w,v,t)\ < Ci\\w\\ H ^ 0t i)\\v\\H^(o,i) (8) 
\A(w,v,t)\ < C 2 || w || H i(o,i)ll«lki(o,i) (9) 

for any t £ [0,T]; 



2. (Garding inequality) there exist constants a > and 7 > such that 

\A(v,v,t)\ > a\\v\\ 2 Hi{01) -7l|w|lia (0il) , (10) 

uniformly with respect to t £ [0,T]. 

Proof. The proof is given in |9J . 

Owing to Lemma [l] one can prove the following assertion [9]. 

Theorem 1. Suppose that u (x) G H^(f2). Then the problem ([7| has an unique 
solution. 

Theorem 2. (Weak maximum principle) Let u 6 L? (0,T; H^(f2)) be such 
that 

1. ff ei 2 ((0,l)x(0,T)); 

2. the inequality 

+ A (u(t),v;t)>0, VveC?(0,l),v>0 

holds for a.e. r £ [0,T]; 

3. u\ t=0 > 0. 

Then u > a.e. in (0,1) x (0, T) . 
Proof. The proof is given in [5J . 

3 Space discretization 

In this section we describe the finite volume approximation of Q. 

Let the interval fl — (0, 1) be subdivided into N intervals Ii := (xi, Xi+i), i = 
0, . . . ,N — 1, with —: x a < x 1 < ■ ■ ■ < x N := 1. For each i = 0, . . . ,N — 1 
we set hi := — Xi and h := maXj=o,...,jv-i hi- We also denote x i+1 / 2 '■= 
[xi + hi/2) for i = 0, . . . , N — 1, £-1/2 : = = 0, Xjv+1/2 : — ^ejv = 1 an d 
J7j := (^-1/2,^+1/2) for i = 0, . . . , N. Finally, we define h := x i+1/2 - ^-1/2 
for i = 0, . . . , N. 

Integrating ^ over the control volumes Qi we obtain JV+1 balance equations 

/ iidx-[x{l- x)p(u)] x *+ 1/2 + [ cudx = 0, i = 0,...,N. 

Multiplying the i-th equation with an arbitrary number Vi and summing up 
the results we get 

N N N . 

^2 uv i dx-^[x{l-x)p{u)] x J+_\'^v i +^ J cu Vl dx = Q. (11) 



For an arbitrary function v £ C(fi) we define the mass- lumping operator 
L : C(Q) -> L^n) by L h v\„ t :=v(xi), i = 0,...,N. 

Therefore, using the operator L^, equation (11) can be written as follows: 



(u(t),L h v) + A h (u, v; t) = 0, Vv £ C{Q), 
A h (w, v; t) := - Yn=o K 1 ~ x )p( w > x , L hv(xi) + (c(t)w, L h v). 

The spatial discretization starts from this equation. Applying the mid-point 
quadrature rule to all terms except the second one we obtain for all v £ C{fl) 



{u{t),v) h - £t [x(l - x)p(u(t),x,t)] x x ^/ / l L h v{xi) + (c(t)u,v) h = 0, 

(W,v)h ■= (L h W,L h v) =J2i=0 W i V i l i, w,v £ C(Q). 

Clearly, we now need to derive approximations of the continuous weighted 
flux density x(l — x) p{u(t) , x , t) , defined above, at the midpoints x i+1 / 2 °f f ne 
intervals Ii for i = 0, . . . , N — 1. 

Case 1 Approximation of p at Xi+1/2 for 1 < i < N — 2. 

Let us consider the following two-point boundary value problem for x £ Ii 

(a l+1 / 2 x(l - x)v' + b i+ i/ 2 v)' = 0; vfa) = Ui,v(x i+ i) = u i+ \, 

where a i+1/2 = a(x l+1/2 ), b l+1 / 2 = b(x i+1/2 ,t). 

Following considerations, similar to those in |6I17| . we obtain 



Pi(u) = b i+1/2 ± X ; +1J / X 'U , i = l,..., N -2, (12) 









(l-x 4+1 ) ( 


Xi 


r 




^ l-Xi 





where ctj = a '^^ and Pi(u) provides an approximation to p{u) at x i+ i. 
Case 2 Approximation of p at x 1 / 2 . 
Now we write the flux in the form 

r \ _ du _ ct 2 

p[u) := ax— h bu, a = a(x) = — (1 — x). 

ax 2 

Note that the analysis in Case 1 can not be applied here because the flux 
degenerates near x = 0. To solve this difficulty, following [2 6 17 20 , we will 
reconsider the flux ODE with an extra degree of freedom in the following form 

(a 1/2 xv' + b 1/2 v)' = C in (0,xi), v(0)=u , v(x x ) = u x , 
where C is an unknown constant to be determined. We obtain 

»(*)-|^(l + ao)) a!3Vl + tlOtaD<0j 



A) 



_ = f S(l + a )x^^ + 6 1/a «o, «o > 0, 



(13) 



Case 3 Approximation of p at Xn-x/2- 
We write the flux in the form 

/ \ = /-, \9u = a 2 

P(U) :=%-l/2 ( 'dx + b N-i/2 u , a [X) = yi- 

The situation is symmetric to Case 2. We can not apply the arguments in Case 1 
to the approximation of the weighted flux density on In-i = (^w-ij #jv) because 
equation Q degenerates at x = xn = 1. However the considerations, given in 
Case 2, should be modified in order to formulate an appropriate two-point BVP. 
Again, we consider the flux ODE with an extra degree of freedom in the following 
form (recall a N -i = b N _ 1 / 2 /a N _ 1 / 2 ) 



((1 - x)v' + au-iv)' = C , x e In-i, 
v(x N -i) = ujv-i, v(x N ) = u N , 



(14) 
(15) 



where Co is an unknown constant to be determined. Integration of ( 14 1 yields 
the first-order linear equation 



(1 — x)v' + cxn-iV = Cqx + Ci, x € In-i, 



(16) 



where C\ denotes an additive constant. We can easily obtain ckn-iUn = Co + Ci 
and this would be legal if we solve the problem with Dirichlet boundary condition 
at x = 1. However, this is not the case and alternatively we multiply (16 1 by 

(1 - i)-""- 1 - 1 



((l-x)- a »-iv) =C (l-x) 
Case 3.1 a^-i > 0, a^-i ^ 1. 



-QJV-I — 1 



x + C x (l-x) 



-Oit-l-l 



(17) 



Integrating (17 1 from xn-i to x <E In-i results in 



(1 - x)- aN - 1 v{x) - (1 - x N ~i)~ aN - 1 v(x N -i) 



Co 



(i-s;r°"-i +1 

— ajv-i+1 



(Co + Ci 



Multiplying both sides of the equation by (1 — x)"™- 1 



(l- a:j v-i)^°"-i + 1 (l-x)'-jv-i 

-OAT-l + l 



-(Co + 61)3^77 + (°o + G i) =^7^ • 

Letting x — > xn — 1 we arrive at v(x^) = Un = ^j"^ 1 and finally 



+ -x)(l- -^=^ 

where u> = Co G R is a free parameter. Therefore 

Pn-i(v) := a N _ 1/2 {x - 1)uj + b N _ 1/2 u N . 

Case 3.2 ajv-i = 1- 

Now we solve the following ODE 



C (l-2;)- 2 a; + C*i(l-x)" 



1 - x, 

Integrating over (xn-i,x),x € In-i, we obtain 
Letting x — > ijv = 1 one gets 



w(xjv) = Mat = C + C\ 



dJV-l 



(1 — x) 1 — X 

v(x) = — r(UjV_i — Mat) + u N + cu(l — x) In 



(1-Xjv_i) ' " l-ZjV-i 

Since 1 — x > we can conclude that this is the result of the limiting process 



Q!JV-1 - > Ij performed on (181. The flux in both cases 3.1 and 3.2 can be written 
in the form 

PN-i(v) = Ojv-i/aOc - l)w + b N -i/ 2 u N . 
Case 3.3 ajv— l = 0. 

Integrating over (xjf-i, x), x £ In-i the following ODE 
v' = C (l - x)- 1 x + Ci(l - x)- 1 

we arrive at 



v(x) = v(x N -i) - C a (x - x N -i) - (C + Ci)(ln(l - x) - ln(l - Sjv-i))- 

The function v(x) is bounded for a; — > xn when Co + Cl = 0. Therefore 
C = T' 1 """ and 



v(x) = U N + (u N -l - U N ) 
The flux has the following form 



1 — x 

1 - X N - 



Pn-i{v) = ajv-iml - ckjv— _ — r H "N-1/2UN, 

n.N-1 

where we used that cxn-i — frjv-i = 0. 
Case 3.4 aw-i < 0. 

This time we integrate from x to xn — 1 and obtain 

= go T^— - + ( C ° + Ci) — 

1 — GLN — 1 CHjV-1 

and using the boundary conditions 



un-i — un n . _ mat-i — UJV N 

C = — (1 - on-i), Gi = ocn —iUn (1 - ajv-ij. 

I-xjv-i l-a:jv-i 

Therefore 

u(as) = ( - — — — J (un-i - un) + u N , 
\l - x N+1 J 

Pn-i(v) = a N -x/i{l - ajv-i)(l - x) h b N _ 1/2 u N 

ftjV-l 

and these are exactly the same results as in Case 3.3. Finally, a reasonable choice 
of the free parameter ui is and 



v(x) 



un + i- x *_ t (un-i - u N ), a N -i<0, 



UN + (i-x^I^n-i (UN-1 - Un), OtN-1 > 0, 



<N-1 



, n I a iv-i/2(l - "jv-i)(1 - a:) , lM , h b N _ 1/2 u N , ajv-i < 0, 

Piv-i(uj = < , llN - 1 

{ ON-1/2UN, ajv-i > 0. 

Let us introduce the finite element space Vh by specifying it's basis {<fii}fLo- 
Following |6I17| we introduce the functions 



— 1 , X t [Xi — i^Xiji 



4>i(x) 



X £ (Xi,X i+ l). 



On the intervals (0,a;i) and (xn-i, 1) we define the linear functions 



Yvy ' \ 0, otherwise; VN KJ | 0, otherwise. 

Next we define the linear functions 4>i(x) and 4>n-i(x) on the intervals (0, x^) 
and (xjy_ 2 , 1) 

l-f- ± ,xe (0,x0, 

(jl-i^-^l-i)" 1 ' ^^(xi,^); 0, otherwise; 



(^^-i) QJV " 2 -(i-i)° N - 2 



(1 — x)/(l — xjv-i), a; € (aijv-i, 1); 0, otherwise. 



Then, for any U/j e V/j, we have the representation Vh(x) — Yli=o v hi fyii 
where Vhi := Vh(xi). Associated with Vh, we introduce the natural interpolation 
operator Ih : G{Q) — > Vh by 

IhVh{xi) := «ft(a;i) = i = 0,...,N. 

Furthermore, using the flux approximations, obtained in Cases 1, 2 and 3 
respectively, we define by ph(u) an approximation to p(u)\j i := Pi(u), i = 
0, . . . , N— 1. Coming back to this motivates the following semi-discretization 
of ^ in the space Vh,: 

(uh(t),v h )h + A h (u h (t),v h ;t) = Vvh G V^, 

where, using the long notation, i.e. p h (v h ,x,t) = p h {vh) with ph(vh,Xi) = pi{v h ) 
for x € Ii 

N 

A h (w h ,v h ;t) := [x(l - x)p h (w h , x, L h v h (xi) + (c(t)w h , v h ) h . 

i=0 

As usual, from |7]) an equivalent ODEs system is obtained by setting suc- 
cessfully Vh = <pi, i = 0, . . . , N: 

U hl {t)k - X l+1/2 (l - X i+1/2 )ph(u h (t), X i+1 / 2 , t) + Xi_i/ 2 ph{Uh{t), 

+c l (t)u h i(t)k = 0, Ci(t) := c(x u t). 

The complete set of equations forms an (N + 1) x (N + 1) system of linear 
ODEs w.r.t. u h (t) := (u h0 (t), u hN (t)) T : 

-M h u h (t) + A h (t)u h (t) = 0, 

where 

M h ■= (Oj,<^)h) i)j=0 = diag(l ,...,l N ), 

A h (t) := (A(^,&;t))J =0 = K-W)5=o- 



4 Full discretization 



Let =: to < t\ < . . . tjy t ■= T be a subdivision of the time interval [0, T] with 
the step sizes At m :— t m+ i — t m > 0, m G {0, . . . , N t — 1}. The fully discrete 
method with parameter 8 G [0,1] for ^ reads as follows: 

Find a sequence U , . . . , U * € V/, sitc/i i/iai for m £ {0, . . . , iVj — 1} 

+^(^ m+1 + (i - ^)c/ m ,^; We) = \/v h e v hl 

where t m+ g :— 8t m+ i + (l — 9)t m = t rrl + 9At m anduoh G Vh is an approximation 
to Uq. By representing the elements U m in terms of the basis {^i}^ 1 of Vh and 
choosing Vh = <t>j 1 j = 0, . . . , N we get the algebraic form 



(i\/r n' m + 1 — i\/r n m \ 
At m ) + 9A ™ +6 < +1 + 0- ~ d )K +6 < = 0, AT := A h (t m ). 

(18) 

The initial condition u° is obtained from the representation of uq^ by means 
of the basis of Vh ■ 

We will show, Theorem jij that the system matrix E/j = {eijjfj^g can be 
reduced to an M-matrix by excluding the first two and the last two equations in 



( 18 ). Therefore, the above problem ( 18 1 is uniquely solvable and our method pre- 
serves the positivity, Theorem [2] (maximum principle), of the numerical solution 
in time. Let us introduce the notations 

oti ( X i \ AOti 1 

\1-Xij - <Pi 

at±i/2 ■= a{Xi±i/2), o i± i/2 : = o(x i±1 / 2 ,t ), c i±1/2 := c(x i±1 / 2 ,t ) 

and write down the elements of the system matrix: 
for i = if ao < then 

eo,o - + 0x 1/2 (l - x 1/2 )(-b™+ e ) + 9l c™+ e , e ,i = 



- (1 - 1?) (x 1/2 (l - ^{-K^) + 10C" 



m+(l-8) 



and if «o ^ then 

eo,o = + tei/ 2 (l - x 1/2 )0.5(a 1/2 - 6™+ 9 ) + WoC", 
co,i = -^1/2(1 - Xi/2)0.5(S 1/2 + b™+ e ), 

+uft(l - 0) (x 1/2 {l - x 1/2 )0.5(a 1/2 + b™+ {1 - 6) )) ■ 



for i — 1 if a < then 

ei,i - + ^+i/ 2 (l - z i+1/2 )&™+Vr A< + e hc? +e , 

ei,o =-0x 1/2 (l-x 1/2 )(-fo™+ e ), 

ei, 2 - -fe i+1/2 (l - Zi+i/a^i^tf+iA"'. 

Fi = u " (zt - (1 - 0) K^a - ^i /a )C/2 _9) ^ a 04 + «icr +(1 - 9) )) 

^1/2(1 -a:i/2)(-6^/ 1 2" 0) ) 
+< 2 (1 - 6)x i+1/2 (l ^i^^J^'Ci^', z = 1 if a < 0; 

and if a > then 

= + ^+1/2(1 - ii+i/2)^/2€ 4 r 

+&r 1/2 (l - x 1/2 )0.5(a 1/2 + fe™+« 2 ) + 0/ lC ™+ e , 
ei,o = -0x 1/2 (l - ar 1/2 )0.5(S 1/2 - 6™+* 2 ), 
ei, 2 - -fe i+ i/ 2 (l - Zi+i/2)^/Vf+iA% 

+*i /2 (i - ^i/ 2 )°- 5 («i/ 2 + + «icr +(1 " e) )) 

+< (1 - 0)z 1/2 (l - x 1/2 )0.5(a 1/2 b^%- 6) ) 

+< 2 (i - 0)z i+1/2 (i - a ;i+1/2 )6^ 1 2 " 0) CiA i ; 

for i = 2,..., TV- 2 



e M = 7^ + ^+1/2(1 - ^+1/2)6^-1/2^' A"' 

+to i _ 1/2 (l - ^_ 1/2 )6™+; 2 ^-M«!r + 9kcT +{1 - 6 \ 

e»,i-i = -^-1/2(1 - ^i-i/2)&^/2 < ^?-i 1 A-i 1 ' 
e M+i = -^+1/2(1 - a; 4+1/2 )6™ H i ;^ 2 ^ 1 Z\ 4 Q % 

*i - «£ - (i - <?) (^+1/2(1 - ^+i/2)Cv2" 9 Vr A' 

+^ V2 (i-^ V2 )ct5r e vr- i A-:i i +^r +(i - e) )) 

+«JS_ 1 (l - 0)x,_ 1/2 (l - ^-i/aJCv^Cl 1 A-I 1 

- 0)^+1/2(1 - x i+1/2 )b^- 9 U?U A ? t ' i = 2, . . . , JV - 2; 



for i — N - 1 if aAr_i > 



e N -i,N-i = ^ + ^-1/2(1 - Xi-i/zW+u^?- 1 + 01n-iC%± 6 v 

ejV-l,JV-2 = -^AT-l/2(l - X N _ 1/2)^1 1/2 j 

eiv-i ;J v - -toi-i/ 2 (l - ^_ 1/2 )6™+; 2 dM^ 1 , 

= «s v _ 1 (A. _ (i _ 0) ( x ._ 1/2( i _ ^.v^Ct^'V?'- 1 ^! 1 

+a;jv-i/2(l - ^Af-i/2)^^ 1 1 / 2 6) + for-ic^ 1- ^)) 
+«hW-2( 1 - 0)2^-1/2(1 - a;jv-i/2)^ 1 1 / 2 fl) 
+n^(l - e) Xi . 1/2 (l ^iCi/^CA-n 

and if ajv-i < then 

+to;v-i/2(l - ^-i/ 2 )0.5(a™ t/2 - &™V2 + OIn-iC^I 
e N -i,N = -Ox N _ 1/2 (l - x N _ 1/2 )0.5(a i _ 1/2 + b™J 1/2 , 
e N -i,N-2 = -Ox i _ 1/2 (l - a) i _i/ 2 )i>™" 1 y 2 ^ a !i 1 4 a li 1 , 

Fn-i = <W-i -(1-9) ( Xi _ 1/2 (l x^blX^U?- 1 ^ 1 

ft \ n c/ = m +( 1 - ) ,m+(l-0) . , m+(l-0)\\ 

+XN-l/2(t - XN-l/2)0-O{ a i-l/ 2 - & i_l/2 +'iV-lC Ar _i I I 

+«S V _ 2 (1 - (9)^-1/2(1 - «i-l/2)Cv 1 2" fl) d 1 ^I 1 

+ u wv( 1 - 0)^-1/2(1 - ^-1/2)0.5(0(^-1/2) + b(x N -i/ 2 )); 
for z = N if Q!Ar_i > then 

&n,n = — h fejv-i/2(l — a; at- 1/2)^-1/2 + 01jmcJ^ +6 ', ejv,iv-i = 0, 



and if on-i < then 
In 



At r 



+ 9x N _ 1/2 (l - x N _ 1/2 )0.5( 



-m+6 



um+9 



■N-l/2 T u N-l/2 



m+6 



-m+0 



sn,n-i — —9x N _ 1 / 2 (l - x N _i/ 2 )0.5(a N _ 1 / 2 - b™_ 1/2 ), 



iTi „,m 

*N - U hN 



( 1 



N 



m+(l- 



(1-9)1 x N _ 1/2 (l - 3^-1/2)0.5(0^1^2 ' + b™ti/ 2 a> ) 



, m +(l-e) 

+ u WV-l( 1 - 6l )a ; w-i/2(l-a;7V-i/2)0.5(a Ar _ 1/2 -b N _\ /2 '). 

Theorem 3. For any given m — 1,2, . . . , N t , if At m is sufficiently small, the 
system matrix of (18 1, Eh, is an M -matrix. 



Proof. Let us write down the scalar form of (18 1: 



B u™ +1 + C u™ +1 = F Q 
A 1 <+ 1 +5 1 <+ 1 +C'i< 2 + 1 = F 1 
A 2 u™ +1 +B 2 u™ +1 +C 2 uZ l 3 +1 =F 2 



where 



^1 = -6>ei, , #1 = 3*^ + 0e ltl , d = -6e h2 , 
Ai = 5; = ^ + 0e M , C, = -0e M+ i, i = 2, 3, . . . , N - 1, 

v4at = —9ejy,N-l, Bjy = 2 At^. @ eN , N > 

Fo = (^t ~ C 1 - ) e o,o) <o + (1 - »)eo,i«M. 
Fi = (1 - 0)ei,otCo + (at " (1 " *Ki) *Ci + (1 - 9)ei,2«JS. 
F, = (1 - ejei.i-xttS-i + (at - (1 - *K*) «JS + (1 " ^)ei,i + i< + i, 
Fy = (1 - 0)e*,JV-i«Jk_i + (53^ - (1 ~ 0)^) 

Let us first investigate the off-diagonal entries of the system matrix A4 = 
—9ei.i-i and C, = — 0ej.j+i. From the formulas for from the above we have 



Gi j > 0, i, j = 1, 2, . . . , N — 1, i ^ j. That is because 



h+1/2 



-l/2<*i 



a i+ i/2Y^T > 0, < n 



Xi+1 1—Xi 



< 1. 



for each i = 1,2, . . . , N — 1 and each 6^+1/2 7^ 0. We have used that 1 — x c * i 
has just the sign of a* and this is also true for b i+ i/ 2 ~ > 0. Now it is clear that 
Ai = —0eij_i and C- L = —Oe^^x are negative. 

We should also note that Bi is always positive since At m is small. The situa- 
tion is different for Bq, Co, Ai, Bi, C\ and An—i, Bn—i, Cv— 1, An, Bn- From 
the first three equations we find 



"h0 — B (> Bo u hl > 



Fa 
Bo 

A = Bi 



m+l 
l hl 



Ai _ Ql.m+1 

A A u h2 > 



^Co, A, 



Fx 



-F , 



B 2 u™ +1 + C 2 u 



h3 



Bo = Bo — 



A 2 C 
A 



L , F 2 = F 2 - ^A 2 . 



It is easily to see that when A > and A = O ( -A— ) then B 2 = O 




for small At m . Therefore B 2 > O (a^J and B 2 > \C 2 \. 

In a similar way one can eliminate u™ N ~^_ 1 and uj^v • a resu lt we ob- 
tain a system of linear algebraic equations with unknowns u™ 2 +1 , . . . , u™jv-2 anc ^ 
coefficients matrix that is an Af -matrix. 

While F 3 , ...,Fn-3 are non-negative, we have to prove if F 2 and Fn- 2 are 
also non-negative. From the formulae for F 2 it follows that when At m is small 
F 2 is non-negative since F 2 = O (ai - ) an d A, Ai are of the same order with 

respect to At m . Fn- 2 is being handled the same way as F 2 and also considered 
non-negative. 

Since the load vector (F 2 , F3, . . . , Fn-3, -Fiv-2) is non-negative and the cor- 
responding matrix is an M-matrix we can conclude that u™ 2 +l , . . . , u™^ 2 are 
non-negative. Finally, using the formulas for u™^ 1 , u™jv_i: u™^ 1 one can 

easily check that they are non-negative too if At m is small. □ 



Remark 1. Theorem [3] shows that the fully-discretized system (18 1 satisfies the 
discrete maximum principle and therefore the above discretization is monotone. 
This guarantees the following: for non-negative initial function uq{x) the numer- 
ical solution u™, obtained via this method, is also non-negative as expected, 
because the price of the option is a positive number. 



5 Numerical experiments 



Numerical experiments, presented in this section, illustrate the properties of 
the constructed method. We solve numerically various European Test Problems 
(TP) with different final (initial) conditions and different choices of parameters. 

1. (TP1). Call option with final condition Parameters: S max = 700, T = 1, 
r = 0.1, a = 0.3, d = 0.04 and E = 400. 

2. (TP2). Call option with cash-or-nothing payoff V(S,T) = H(S - E), S £ 
(0, oo), where H denotes the Heaviside function. Parameters: S max — 700, 
T=l,r = 0.1, cr = 0.3, d = 0.04 and E = 400. 

3. (TP3). Call option with final condition Parameters: S max = 700, T = 1, 
r = 0.1 + 0.02sm(10Ti), d = 0.06a; and E = 400. 

4. (TP A). Portfolio of options. Combinations of different options have step final 
conditions such as the 'bullish vertical spread' payoff, defined in (|3j. In this 
example, we assume that the final condition is a 'butterfly spread' delta 
function, defined by 

l, Se(S u S 2 ), 
V(S,T) = { -l,5e (S 2 ,S 3 ), 
0, otherwise, 

which arises from a portfolio of three types of options with different exercise 
prices [!§] . 

In the tables below are presented the computed C and L 2 mesh norms of the 
error E — Uf, . — u by the formulas 



|£|| c =max||<*-^||, ||£|| i2 



N „ 



\ ^ 

\ i=0 

while the rate of convergence (RC) is calculated using double mesh principle 

RC = log 2 (E N /E 2N ), E N = \\u%-u N \\, 

where ||.|| is the mesh C-norm or _L2-norm, u N and are respectively the exact 
solution and the numerical solution computed at the mesh with N subintervals. 

In Table [T] we show the convergence and the accuracy of the constructed 
scheme for 9 = 0.5, the weight parameter with respect to the time variable, 
when we numerically solve the model problem with the known exact solution 
u(x,t) = exp(x — t) and initial data uq(x) = exp(x). We choose this function 
because it's feature is similar to that of the exact solution to the call option 
problem. The results, corresponding to problems TPl and TPS with At m = 
0.001, m = 0, . . . , N t - 1, axe listed in Table [l] 

In Table |2]the results are obtained by computations on a power-graded mesh 
for the same values of the parameters and exact solution. This mesh takes into 



Table 1. 



N 




TP1 






TPS 




tttJV 


RC 


E? 


RC 




RC 


E» 


RC 


80 


3.4547e-3 




2.8008e-4 




4.8049e-3 




3.9144e-4 




160 


1.7292e-3 


0.998 


9.9136e-5 


1.498 


2.4046-3 


0.998 


1.3853e-4 


1.498 


320 


8.6503e-4 


0.999 


3.5068e-5 


1.499 


1.2028e-3 


0.999 


4.9000e-5 


1.499 


640 


4.3256e-4 


0.999 


1.2400e-5 


1.499 


6.0147e-4 


0.999 


1.7326e-5 


1.499 



account the degeneration at the both ends of the interval (0, 1) and is given by 
(in the current case p = 2) 

z i+1 =Xi+(i (2E!I / o 2 ^ 1/P ) ! ', i = l,..,N/2 

Xi+i = a* + ((N + 1 - i) (2 Y,^lo i p )~ 1/P y ,i = N/2 + 1,.., N. 

The time step At m is chosen such that At m = mino<i<Ar h(i),rn — 0, . . . , N t — 1 
with T = 0.1. 



Table 2. 



N 




TP1 






TPS 






RC 


E? 


RC 




RC 


E$ 


RC 


20 


7.1539e-4 




3.6478e-4 




6.2631e-4 




3.9140e-4 




40 


1.8803e-4 


1.927 


9.5247e-5 


1.947 


1.6496e-4 


1.924 


8.3410e-5 


2.230 


80 


4.8181e-5 


1.964 


2.4372e-5 


1.966 


4.2261e-5 


1.964 


2.1344e-5 


1.966 


160 


1.2195e-5 


1.982 


6.1671e-6 


1.982 


1.9695e-5 


1.982 


5.4011e-6 


1.982 



In the last table we compute the solutions of the original models TP2 and 
TP3. As exact solution we use numerical solution on a very fine grid, i.e. N = 
5120 with At m = 0.0001, m = 0, . . . , N t - 1. The results are given in Table [3] 

The test parameters as well as the obtained results (see also Fig(T|4]) are 
similar to those in [20] as expected. 

6 Conclusion 

In this article we present a fitted FVM for the generalized Black-Scholes equation 
(JTJ. The method is applicable to more general Black-Scholes models, for example 
when a = u[S, t) and r — r(S, t). We may also use any interval (0,1), I > (here 



Table 3. 



N 




TP2 






TPS 




rpN 


RC 


E? 


RC 


pJV 
E x 


RC 


E? 


RC 


80 


2.9139e-7 




1.1120e-7 




2.6809e-3 




2.1714e-4 




160 


9.9137e-8 


1.555 


2.8407e-8 


1.968 


1.3205e-3 


1.021 


7.4759e-5 


1.538 


320 


5.0468e-8 


0.974 


7.3863e-9 


1.943 


6.3931e-4 


1.046 


2.5443e-5 


1.554 


640 


2.5449e-8 


0.987 


1.9726e-9 


1.904 


2.9842e-4 


1.099 


8.3735e-6 


1.603 


1280 


1.2689e-8 


1.003 


5.418e-10 


1.864 


1.2791e-4 


1.222 


2.5343e-6 


1.724 




Fig. 1. Test Problem 1 Fig. 2. Test Problem 2 



we took I — 1 for simplicity) to solve the transformed problem. The main ad- 
vantage of the developed numerical algorithm is reduction of the computational 
costs as well as positivity-preserving. 

In a forthcoming paper we study the stability and the convergence of the 
proposed finite volume method. 
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